Multiple light scattering in laser particle sizing 


Alexander A. Kokhanovsky and Reiner Weichert 


Multiple light scattering is an important issue in modern laser diffraction spectrometry. Most laser 
particle sizers do not account for multiple light scattering in a disperse medium under investigation. 
This causes an underestimation of the particle sizes in the case of high concentrations of scatterers. The 
retrieval accuracy is improved if the measured data are processed with multiple-scattering algorithms 
that treat multiple light scattering in a disperse medium. We evaluate the influence of multiple light 


scattering on light transmitted by scattering layers. 


The relationships among different theories to 


account for multiple light scattering in laser particle sizing are considered. © 2001 Optical Society of 


America 
OCIS codes: 


1. Introduction 


Small-angle laser diffraction is often used for deter- 
mining particle size distributions (PSD’s) in disperse 
media.! The retrieval schemes in this case are 
based on the strong dependence of Fraunhofer dif- 
fraction patterns of large particles on their sizes. In- 
deed it follows for the scattered light intensity in the 
case of monodisperse ensembles of large spheres? 
that 


J,°(a0)a” 


I(0,a)=N OR 


IV, (1) 


where J, is the Bessel function; a is the radius of 
particles; a = ka, where k = 27/\, where \ is the 
wavelength of incident light; 6 is the scattering angle, 
I, is the intensity of the incident light, N is the num- 
ber of particles per unit volume, V is the volume of a 
scattering medium, and R is the distance to the ob- 
servation point. The intensity J(6, a) is proportional 
to the square of the Bessel function and therefore has 
several minima and maxima. The positions of min- 
ima and maxima depend on size parameter a. The 
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first minimum occurs at scattering angle 0 = A/a, 
where A ~ 3.832. 

We can account for the polydispersity of a disperse 
medium by averaging the scattered light intensity 
[see Eq. (1)] with the PSD qola): 


as | Maide (2) 


0 


Equation (1) is valid under the following assump- 
tions: 


a > 1, 2a\m - 1| > 1, (3) 
a = 1 (4) 
ip ? 
4 > 1 (5) 
d ? 


where m = n — ik is the refractive index of the 
particles, L is the geometric thickness of a scattering 
layer, / is the mean photon free path length, d is the 
diameter of the particles, and D is the distance be- 
tween particles. Note that, owing to assumption (3), 
the Fraunhofer diffraction pattern of a particle coin- 
cides with the Fraunhofer diffraction pattern for a 
black screen with radiusa. The task of the retrieval 
algorithm of the laser diffraction spectrometry is to 
find the radius (or radii distributions) of such screens. 
Equations (4) and (5) state, respectively, that the 
multiple scattering? and correlation effects*> are neg- 
ligible. 

In many problems, such as those encountered in 
particle sizing for Diesel injectors and heavy fuel oil 
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atomizers, the correlation effects can be neglected (D 
>> d), but condition (4) is often not fulfilled. Thus 
one needs to account for multiple light scattering and 
obtain a more general formula than the simple single- 
scattering approximation presented by Eq. (1). 

This could be done by different means. Interest- 
ing approaches to dealing with the problem have 
been proposed by Felton et al.6 and Hirleman.”® Op- 
tical particle-sizing techniques based on the use of 
the radiative-transfer equation? have been studied by 
Belov et al.,2 Vagin and Veretennikov,!° and 
Schnablegger and Glatter.1 

We show that analytical solutions’!! are identical 
at small angles and can be mutually derived from one 
another. The simplified equation for the calculation 
of small-angle light-scattering patterns of multiple 
light-scattering media with large spherical particles 
is derived. 


2. Multiple Light Scattering 


A. Different Forms of Solutions 


The radiative-transfer equation (RTE) for a plane- 
parallel light-scattering layer can be written in the 
following form?: 


dI (r, 3) 


ù 
cos dz 


= =17; ù) 


+ > f I(r, 9’)p(@)sin 9'd9’, (6) 
0 


where 9 is the observation angle; wo = Osca/Text is the 
single-scattering albedo, where o,,, is the scattering 
coefficient and o,,; is the extinction coefficient; p(6) is 
the azimuthally averaged phase function; T = o,,,L is 
the optical thickness, where L is the geometric thick- 
ness of a layer; and IJ is the light intensity. The 
scattering angle 0 is related to angles ð, 3’ by the 
following simple equation?: 


6 = arccos[cos ð cos ð’ + sin ð sin ð’ cos (ð — 8’) ] 
(7) 


where ò and 0’ are azimuths. 

The specific form of the RTE presented by Eq. (6) is 
valid only at the normal illumination of a light- 
scattering plane-parallel layer bya plane wave. The 
transmitted diffuse light intensity does not depend on 
the azimuth in this case owing to the symmetry of the 
problem. 

We consider the angular dependence of function 
I(t, Ò) at small values of the observation angle ð. 
Thus it is assumed that 


a Ma, 9) _ all, 8) 
a dt dtr’ 


(8) 


where we accounted for the approximate equality: 
cos ð ~ 1, which is valid as 8 —> 0. This approxi- 
mation is called the small-angle approximation? of 
the radiative-transfer theory. Equation (6) under 
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approximation (8) can be solved analytically. The 
answer is well known (see, e.g., Refs. 3 and 12): 


I feo) 
Ir, 8) =—> E(t) P,(cos 9), (9) 
2T 7 
where 
= sg al 
Ei) = | + 5) expl-( — oop), (10) 


1 (7 
p= A p(0)P;(cos @)sin 6dé. (11) 
0 


Here P, (cos 0) is the Legendre polynomial. 

Our prime interest here is the incoherent intensity 
of multiply scattered light. Thus we remove the co- 
herent component? (or direct light), 


I 
I, exp(—t) 8(cos ð — 1) = a exp(—7) 


T 


= 1 
x > | + z] os ð), 
j=0 2 


(12) 
from Eq. (9) and obtain 
Io = — 
I'a, 8) =— > Z, 0)P;(cos 9), (13) 
2r j=0 


where 


1 
=; 0) = f + 3) e “[exp(wop,t) — 1]. (14) 


Equations (13) and (14) were obtained by Hartel!” 
and used in Ref. 11 for solution of the inverse prob- 
lem. Note that numerical calculations of coefficients 
p; [see Eq. (11)] are complicated in the case of large 
particles because the phase function is strongly 
peaked in the forward direction (0 = 0) ata > À. 


Let us derive the integral form of Eq. (13). It fol- 
lows as 0 — 0: 
ee | 
P,(cos 0) > Jo (i + 9) ; (15) 


Thus we can introduce the following function, which 
follows from Eqs. (11) and (15): 


1 
so =| 


0 


o 


 p(©)J(00)0d9, (16) 


where o =j + ¥% and we extend the upper limit of the 
integration in Eq. (11) to infinity. This is possible 
owing to a sharp peak in the phase function p(0) [see 
Eq. (1)] at scattering angle 0 ~ 0. Also we use the 
continuous function g(c) instead of discrete numbers 


p; [see Eq. (11)]. It follows from Eq. (13) and the 
Euler sum formula 


Sri) [ roa 


that 


T' (t, 9) = a e” T {exp[Tw)g(c)| — 1}Jo(od)odo. 
0 
(17) 


Equation (17) is more convenient for applications 
than Eq. (13). Note that the simple form of Eq. (1) 
allows for the analytical integration in Eq. (16) (see 
Subsection 3.A). 

Equation (17) was used in Refs. 9 and 10 for solu- 
tion of the inverse problem taking into account mul- 
tiple light scattering. Thus it is obvious that 
approaches to optical particle sizing in Refs. 9, 10, 
and 11 are based on asymptotically equivalent for- 
mulas [see Eqs. (13) and (17)]. Note that Eq. (17) 
can be obtained directly from Eq. (6) by the use of 
Fourier transform techniques.!* 

Let us now consider the expansion of the exponent 
in Eq. (17): 


œ 


T'w'8"(0) 


exp[To)g(0)] = >) —— (18) 
n=0 n: 
It follows from Eqs. (17) and (18) that 
I o0 n n 
T'i, 8) =e" S E 00), (19) 
Qa n=1 n! 
where 
H,(8) = l "(oJ (oð)odo. (20) 
0 


Note that the successive scattering order expansion 
[Eq. (19)] coincides14 with Hirleman’s solution [see 
Eq. (33) in Ref. 7] and provides a new way to calculate 
the function H,,(9). 


B. Single-Scattering Approximation 


It follows from Eq. (19) at n = 1 and small values of 
t (the single-scattering approximation) that 


I œ 
T' (t, 9) = a wo | g(o)Jo(od)odo, (21) 
T 
0 
or using the equality Two = Oscal yields 
Io 
I(r, Ò) Se Oscalp (9), (22) 
20 
where [see Eq. (16)] 
p(s) = | 2g(o)Jo(od)odo (23) 
0 


is the phase function. On the other hand, it is well 
known? that the phase function of the Fraunhofer 
diffraction process has the following form: 


4.J ,°(Sp) 
p(s) = (24) 
Ù 
Thus for the intensity within the framework of the 
single-scattering approximation one can obtain 


I(t, 9) = p(®)a2NLIo, (25) 


where we use the equality oga = 27Na”, which is 
valid at a > i.3 

It follows that Eq. (25) [see also Eq. (24)] coincides 
with Eq. (1) atL = V/(4R*). Different multipliers [L 
and V/(4R°)] are due to the different symmetry of the 
problem associated with Eqs. (1) and (6). Note that 
normalized intensities i(t, ©) = I(t, ©)/I(r, 0) defined 
by Eqs. (1) and (25) coincide. 

Thus one can see that Eqs. (13), (17), and (19) are 
special cases of the solution of the RTE [Eq. (6)]. 
They have already been used for development of dif- 
ferent inversion schemes within the framework of 
laser diffraction spectrometry.®:9:11 


3. Multiple Fraunhofer Diffraction 


A. Monodisperse Media 


The main task in this section is to study Eq. (17) in 
more detail for the multiply scattered light intensity 
I(t, ò). This formula can be applied to many kinds of 
disperse media with different types of forward 
peaked phase function, including Rayleigh—Gans 
phase functions.* 

Let us now consider a special case of Eq. (17) with 
the phase function p(6), defined in Eq. (24). It fol- 
lows for the Fourier—Bessel transform of the phase 
function (see Fig. 1) that 


2 
—[arccos(z) — z(1 — 2)"?] z<1 
T 


g(z)= , (26) 


0 g>1 


where z = o/n and y = 2a. The single-scattering 
albedo wọ is equal to 0.5 within the framework of the 
Fraunhofer approximation. Thus it follows from Eq. 
(17) for the intensity of the transmitted diffused light 
that 


Io 
I(r, 8) =e" 
S m | 


0 


“ [exp 2S] = 1edd 
(27) 


or 


I(t, 8) = Ce "7 f Gia = 1 bee, (28) 


where b = 2kað, C = [(2a”)/m]J, and the function 
g(z) is presented by Eq. (26) and Fig. 1. One can 
observe the results of calculations of the normalized 
intensity i(t, Ò) = I(t, ©)/I (7, 0) obtained with Eq. (28) 
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Fig. 1. Fourier—Bessel transform of the Fraunhofer phase func- 
tion [Eq. (24)] obtained with the exact solution [Eq. (26)] and the 
approximate formula [Eq. (34)]. 


for different values of the optical thickness 7 as a 
function of parameter b in Fig. 2. It follows that 
multiple light scattering causes broadening of the 
angular intensity distributions. The same effect 
could be due to the decrease in particle size [see Eq. 
(1)]. Thus multiple-scattering broadening is respon- 
sible for the underestimation of the particle size by 
laser diffraction spectrometers, designed for thin 
scattering layers, if they are used for large values of 
the optical thickness. 

Let us define the dimensionless halfwidth param- 
eter h by 


i(b = h) = 0.5. (29) 
It can be calculated from Eq. (28) for different t. The 
value of h increases with the optical thickness (see 
Fig. 3 and Table 1). Note that Table 1 can be used to 
estimate the radii of monodisperse particles from 


T 
0.01 — 
f cess 


Fig. 2. Dependence of the normalized intensity on parameter b = 
2kaù obtained with Eq. (28) at different values of optical thickness 
t = 0.01, 1, 2, 5, 7. 


1510 APPLIED OPTICS / Vol. 40, No. 9 / 20 March 2001 


4.4 


42h : : J 


half-width 
w 
& 
T 
i 


3.2 i i i i i i 
0 1 2 3 4 5 6 7 
optical thickness 


Fig. 3. Dependence of half-width / of the angular spectrum of the 
transmitted light on optical thickness r. 


measurements of optical thickness 7 and observation 
angle ðo at which i(ðọ) = 0.5: 


h(t) 
a — 
Ato 


N: (30) 


The value of h in Eq. (30) should be selected from 
Table 1, depending on optical thickness. Also one 
can use h = 3.23614 + 0.07687 + 0.009377”, which 
was obtained by fitting the data in Table 1 (see Fig. 
3). The same method can be applied for estimation 
of the mode radius of narrow particle-size distribu- 
tions. However, note that the values in Table 1 ig- 
nore light rays reflected and transmitted by particles. 
They account only for the diffraction component of 
the total scattered-light intensity. 

It is useful to estimate the upper limit of optical 
thickness t*, where multiple scattering can be ne- 
glected. It follows from Fig. 2 that the influence of 
multiple light scattering is larger for larger values of 


Table 1. Dependence of the Half-Width Parameter h of the Angular 
Distribution of Transmitted Light on Optical Thickness 7 for Disperse 
Media with Monodispersed Spheres 


T h 
0.01 3.23 
0.5 3.28 
1.0 3.32 
15 3.38 
2.0 3.43 
2.5 3.49 
3.0 3.55 
3.5 3.62 
4.0 3.69 
4.5 3.77 
5.0 3.85 
5.5 3.94 
6.0 4.03 
6.5 4.14 


0.1 E 


optical thickness 


0.001 i F a 
0.1 1 10 100 
obscuration, % 


Fig. 4. Dependence of the optical thickness on obscuration. 


parameter b = 2kað. One can introduce the follow- 
ing criteria on the basis of data in Table 1: 


h(r*) 
h(0.01)” er 


where the value of € depends on the tolerable error for 
the retrieval scheme. Clearly, the value of h at t = 
0.01 corresponds to the case of single scattering; e.g., 
if one allows for the difference in half-widths £ = 1%, 
it follows that 7* < 0.5 to fulfill Eq. (81). Malvern 
Instruments advises that the upper limit of the ob- 
scuration p = [1 — exp(—7)]| X 100% should be 50% to 
avoid multiple light scattering.!5 The optical thick- 
ness T* is equal to 0.7 in this case (see Fig. 4) and € ~ 
2% [see inequality (31) and Table 1]. Bayvel et 
al.16.17 found that it is possible to obtain reliable re- 
sults of the inversion even at u = 90% (or at T* = 2.3, 
see Fig. 4). One can find that the value of € is equal 
to 7% in this case. 

In conclusion, Eq. (28) [see Eq. (26) as well] pro- 
vides a simple and useful tool for studies of the con- 
tribution of multiple light scattering to Fraunhofer 
diffraction patterns of scattering media with mono- 
disperse spheres. 


B. Polydispersions 


It is important to generalize Eq. (27) for the case of 
polydispersions. It could be done by averaging opti- 
cal thickness t and the Fourier—Bessel transform of 
the phase function g(z) [see Eqs. (1) and (26) | with the 
PSD q,(a) (Ref. 20): 


ge | A 82) 
" g(o/2ka)a?qo(a)da 
(g(o)) = =— I (33) 


| | a’qo(a)da 


where B = o/2k and t(a) = 2ma?NL (Ref. 3) for large 
particles. One should substitute the values of Eqs. 
(82) and (88) [instead of t and g(o/m)] in Eq. (27) to 
obtain the expression for the intensity of the light 
transmitted through a polydisperse medium. 

Thus it follows that the only complication in Eq. 
(27) for a polydispersed medium is that function g(c) 
is not described by the simple Eq. (26). One should 
use the ratio of integrals [Eq. (33)] in this more com- 
plex case [see Eq. (32)]. The optical thickness of a 
polydisperse medium is determined by the integral fő 
a’qo(a)da. The integrals in Eqs. (32) and (33) can be 
evaluated analytically for some special cases. 

For example, it could be done by expansion of the 
function g(z) [see Eq. (26)] at small values of z: 


3 
g(z)=1+> aa (34) 
j=0 
where sy = —(4/m), sı = 2/(87), sg = 1/(107), sg = 
1/(287). Here we neglected terms of the order of z® 
and higher. The error in the replacement of Eq. (26) 
by Eq. (34) is small (see Fig. 1). The largest error is 
located atz ~ 1. However, the value of g(z) is small 
at small values of 1 — z[g(1) = 0] [see Eq. (26)], and 
the correspondent contribution to Eq. (27) can be ne- 
glected. 
It follows from Eqs. (33) and (34) that 


(g(o)) =1+ >, sj + 1, B), (35) 


where 


o(n, B) = B" ~~ . (36) 


a’qo(a)da 
0 


The special function y(n, B) in Eq. (36) can be calcu- 
lated analytically for many types of PSD. For exam- 
ple, it follows for the gamma PSD qola) = Aa” 
exp(—ap/d), where A is the normalization constant 
(fò Yo(a)da = 1), ap is the mode radius, and p is the 
half-width parameter: 


(g(o)) = 1+ È sj + 1, P), 


as 
2+p-n 


| y exp(—y)dy 
oln, B) = B” — : (37) 
| y?*" exp(—y)dy 


0 


where B = (Apo)/(471a). 
3 one can obtain 


From Eq. (37) at p >n — 


nT t+ 3 —n)[1 -Plu + 3 — n, B)] 


(38) 
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where 
Ta) = | t exp(—t)dt (39) 
0 
is the gamma function, and 


Plana ho | "e exp(—t)de (40) 
0 


is the incomplete gamma function. It follows at in- 


teger values of u that 


p-lal 
T(w)=(u- 1), — P(u, B) = 1- exp(—B) > ; 
1=0 & 
(41) 
and [see Eq. (38)] 


Tw +3—n) gp 


X = 42 


Tuts) A l 


o(n, B) = B” exp(—B) 


Thus for the function (g(c)) within the framework 
of the approximation being considered one can obtain 


3 p-2j+1 
(elo) =1+ > > OF eae 
j=0 1=0 
x exp(—B) hea?) (43) 


Tíu + 3rd + 1)" 


Note that optical thickness (t) [Eq. (32)] is related 
to the volumetric concentration of particles c by a 
simple formula: 


1.5cL 
(1) = ; (44) 
Aeff 
where 
| a°qo(a)da 
a= — (45) 


| a’qo(a)da 


is the effective (or Sauter) radius and c is the fraction 


of a unit volume, occupied by particles. It follows 
from Eq. (45) for the gamma PSD that 
3 
Aeff = Ao [ =F | . (46) 
u 


Equation (44) can be used for estimation of the con- 
centration of particles c from measured values of L, 7, 


Qef- 
Thus it follows from Eqs. (27) and (44) [see Eq. (43) 
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normalized intensity 


TESS oe 
Cee iene 
ee EENE ee 
<4 sae eS haere ae ol 
3.5 4 45 5 


0 L L L 
0 0.5 1 15 


2 25 3 
observation angle,degrees 
Fig. 5. Dependence of the normalized intensity on the observa- 
tion angle obtained with the Mie theory (Mie) and Eqs. (47) and 
(43) for the gamma PSD f(a) = Aa® exp(—1.5a) at wavelength \ = 
0.628 um, the geometric thickness of a scattering layer L = 4 cm, 
and different volumetric concentrations of particles c = 0.000001, 
c = 0.00001, c = 0.0001, c = 0.0002, c = 0.0005, c = 0.0007. 


as well] for the intensity of the diffused transmitted 
light in the small-angle range that 


l Io 3cL \ f” 3cLig(c)) 
ee 20 ex on) J [ex ee 7 | 


0 


x JIo(od)odo. (47) 


Results of calculations of normalized intensities 
i(8) = [1(8, 7)]/[Z(0, +)] with Eq. (47) for polydisper- 
sions with the PSD q,(a) = Aa” exp(—1.5a) are pre- 
sented in Fig. 5 for different volumetric 
concentrations of scatterers c and a wavelength of 
0.628 pm. The geometric thickness of a layer and 
the effective radius of particles were equal to 4 cm 
and 6 wm [see Eq. (46)], respectively. 

It follows from Fig. 5 that the increase in concen- 
tration of particles causes a broadening of the angu- 
lar distribution. This is the same effect as for 
monodisperse media. We also present data for the 
normalized intensity calculated with the Mie theory 
for the same PSD, \ = 0.628 pm and n = 1.33 (water 
drops) in Fig. 5. The small difference between the 
Mie curve and our curves at c = 10°, 10° (7 = 0.01, 
0.1) is related to the error of approximations in Eqs. 
(1) and (48). 


5. Conclusion 


The widely accepted approach to optical particle siz- 
ing is based on the single-scattering approximation. 
However, many natural media (e.g., atmospheric 
aerosols and clouds) are optically thick, and one 
should account for the multiple light scattering 
events in the retrieval schemes. This is important 
for laboratory measurements as well, because some 
disperse media (e.g., in the solid state) cannot be 
diluted. 

We have found that different optical particle-sizing 
techniques’-!! have the same ground, namely, the 


small-angle approximate solution [Eq. (17)] of the 
radiative-transfer Eq. (7). Equation (17) can be ap- 
plied for the solution of both direct and inverse prob- 
lems of disperse media optics in multiple-scattering 
conditions. This formula was transformed to Eq. 
(47) for the simplification of calculations of the small 
angle transmitted intensity in the case of light- 
scattering media with polydisperse large spherical 
particles. 
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